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BACKGROUND OF THE INVENTION 

Field of the Invention 

[0001] This invention relates generally to the field of seismic prospecting and, more 
particularly, to a method for deriving reservoir lithology and fluid content by stochastic 
5 inversion of seismic data. 
Background of the art 

[0002] In the oil and gas industry, seismic prospecting techniques are commonly used to 
aid in the search for and evaluation of subterranean hydrocarbon deposits. Generally, a 
seismic energy source is used to generate a seismic signal which propagates into the earth 
10 and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between 
underground formations having different elastic properties). The reflections are recorded 
by seismic detectors located at or near the surface of the earth, in a body of water, or at 
known depths in boreholes, and the resulting seismic data may be processed to yield 
information relating to the geologic structure and properties of the subsurface formations. 

15 

[0003] The goal of seismic data processing is to extract from the data as much 
information as possible regarding the subsurface formations. Data processing techniques 
have been developed which typically permit the geologic structure of the subsurface 
formations to be determined with a great deal of accuracy. However, to date, efforts to 
20 develop techniques for deriving the fluid content of the subsurface formations have met 
with only limited success. 
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[0004] It is well known by persons skilled in the art of seismic prospecting that the 
presence of hydrocarbon accumulations in a subsurface formation can have a significant 
effect on the velocity of propagation of compressional waves (P-waves) through that 
formation. This effect is the basis of the so-called "bright spot" phenomenon in which an 
5 anomalously high reflection amplitude on a seismic section is an indication of the 
presence of hydrocarbon accumulations, particularly natural gas, in the formation. 
Unfortunately, the bright spot phenomenon is susceptible to error because many seismic 
amplitude anomalies are not caused by hydrocarbon accumulations, or they are caused by 
hydrocarbon accumulations which are low in total saturation and often non-commercial. 
10 For this reason, wells drilled on such bright spots often encounter either no reservoir 
sands at all (and, therefore, no hydrocarbons), or if the sands are present, no 
hydrocarbons or only low saturations of hydrocarbons. 

[0005] One technique which may be useful for this purpose is amplitude variation with 
15 offset ("AVO") analysis. In AVO analysis, measurements of P-wave reflection 
amplitudes with different angles of incidence are used to attempt to determine 
compressional wave (P-wave) velocity, shear wave (S-wave) velocity, density, and 
Poisson's ratio for each subsurface layer suspected of containing natural gas. Knowledge 
of these subsurface elastic properties can be used to predict whether or not natural gas 
20 accumulations are present. See e.g., Ostrander, W. J., "Plane-wave reflection coefficients 
for gas sands at non-normal angles of incidence," Geophysics, v. 49, pp. 1637-1648, 
1984, for a discussion of AVO analysis. Ostrander proposes a method for using AVO 
analysis to distinguish between gas-related amplitude anomalies and non-gas-related 
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amplitude anomalies. However, Ostrander admits that distinguishing between low gas 
saturation and full saturation can be very difficult. 

[0006] AVO techniques have been the subject of a number of prior patents. For example, 
5 U.S. Pat. No. 4,858,200 to Goins ("the Goins '200 patent") discloses a method for 
determining the presence of hydrocarbons in subsurface geological formations by 
comparative assessment of P-wave and S-wave reflection data. The S-wave reflection 
data is estimated from the P-wave data using variations in the amplitude of the gathered 
P-wave data with source-receiver offset. Two related patents, U.S. Pat. Nos. 4,858,201 
10 to Goins et al. ("the Goins '201 patent), and 4,858,202 to F/fc/ietal. describe two 
different methods which can be used for obtaining S-wave data from common depth 
point gathered P-wave traces. 

[0007] U.S. Pat. No. 4,817,060 to Smith discloses a process for directly detecting the 
15 presence of hydrocarbons from seismic data. First, the P-wave and S-wave reflectivities 
are extracted from the data on a trace-by-trace basis. The P-wave reflectivity is then 
determined as a function of the S-wave reflectivity and the result is subtracted from the 
extracted P-wave reflectivity to define a fluid factor which is indicative of the presence 
of hydrocarbons. 

20 

[0008] U.S. Pat. No. 5,001,677 to Masters discloses methods for processing and 
displaying seismic data to emphasize potential hydrocarbon bearing strata. These 
methods treat measured attributes from the seismic data as components of a vector, 
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estimate a background vector which represents uninteresting geologic behavior, and then 
form at least one new attribute which quantifies departures from this uninteresting 
behavior. 

5 [0009] The end products of these prior art AVO processes usually are predictions of the 
P-wave and S-wave reflectivities for the target location. Although some of these prior 
processes have recognized the desirability of also determining the density reflectivity 
(see e.g., the patents to Smith and Masters cited above), none has disclosed a method for 
successfully doing so. 

10 

[0010] Another technique which may be useful in discriminating between different 
lithologies and fluid saturations is pre-stack inversion based on either a one-dimensional 
(ID) or two-dimensional (2D) model of the earth's subsurface. See e.g., Symes, W. W, 
and Carazzone, J. J., "Velocity inversion by differential semblance optimization," 
15 Geophysics, v. 56, pp. 654-663, 1991, and Liao, Quingbo. and McMechan, G. A., 

"Multifrequency viscoacoustic modeling and inversion," Geophysics, v. 61, pp. 1 371- 
1378, 1996. 



[0011] As will be well known to persons skilled in the art, seismic inversion is a process 
20 for deriving a model of the earth's subsurface from seismic reflection data. First, the 

process attempts to extract information regarding the elastic properties of the subsurface 
from the data. This information is then used to construct a mathematical or physical 
model of the earth's subsurface, and synthetic seismograms are generated based on the 
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model. If the synthetic seismograms do not compare favorably to the data, appropriate 
adjustments are made to the model, and new synthetic seismograms are generated for 
comparison with the data. This process repeats until the synthetic seismograms 
generated from the model approximate the actual data. The model is then accepted as 
5 accurate. 

[0012] Pre-stack inversion processes typically attempt to estimate both the background 
P-wave velocity model and the contrast in various elastic parameters (P-wave velocity, 
S-wave velocity, and density) and, therefore, are non-linear. Thus, these techniques are 
10 extremely complex. 

[0013] United States Patent 5,583,825 to Carrazzone et al. describes a method for 
deriving reservoir lithology and fluid content for a target location from pre-stack seismic 
reflection data. The method uses inversion of pre-stack seismic reflection data for both 

15 the target location and a calibration location having known subsurface lithology and fluid 
content to derive the subsurface lithology and fluid content at the target location. The 
inversion process is preferably a viscoelastic inversion to account for the effects of 
friction on seismic wave propagation. The results of the inversion process are a set of 
subsurface elastic parameters for both the target and calibration locations. Relative 

20 magnitudes of these subsurface elastic parameters are compared, together with the known 
subsurface lithology and fluid content at the calibration location, to derive the subsurface 
lithology and fluid content at the target location. 



GTI-1539 



[0014] The method of Carrazzone , while giving results superior to those in earlier 
techniques, still requires an inversion of seismic data and still carries out a two step 
procedure. In the first step, an inversion of the pre-stack seismic reflection data is 
carried out to determine the selected set of elastic parameters at each of a plurality of 
5 points in the models of the subsurface target and calibration locations. In the second 
step, the relative magnitudes of the elastic parameters for the subsurface target and 
calibration locations are compared; and using the results of the comparison and the 
known lithology and fluid content at the subsurface calibration location the lithology and 
fluid content at the subsurface target location are derived. 

10 

[0015] A problem with all of the prior art methods is the two step procedure, explicit or 
implicit, used for obtaining fluid properties. There inversion of seismic data to obtain 
reflection coefficients (or elastic parameters) is by itself difficult. The second step of 
determination of fluid properties from reflection coefficients requires an inversion 
15 procedure that is very sensitive to the unknown parameters being determined. A variety 
of parameters must be used and some of these parameters must be obtained outside the 
inversion itself. It would be desirable to have a robust method of determination of fluid 
parameters of subsurface formations that also takes into account the relative uncertainty 
in knowledge of subsurface rock formations. The present invention satisfies this need. 

20 

SUMMARY OF THE INVENTION 
[0016] The present invention is a method for determining a parameter of interest of a 
fluid in a subsurface region of earth formations. Seismic survey data are obtained over 
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the subsurface region and a horizon of interest is identified. One or more seismic 
attributes of the horizon of interest are measured and a first probability density function 
(PDF) for the seismic attribute(s) is determined. A rock properties model is defined 
consistent with expected values and a plurality of perturbations of the model is obtained. 
5 For each of the plurality of models and a trial value of a fluid property, a PDF of the 
seismic attribute(s) is determined. By comparing the first and second PDFs, an estimate 
of the fluid property is made. .When either the number of seismic attributes or the 
number of fluid properties is more than one, the PDF is multivariate Any one of the 
commonly used seismic survey types (P-P, P-S, S-S, and S-P) may be used. 

10 

[0017] Any of the commonly used seismic attributes may be measured. These include 
impedance, amplitude, a reflectivity, density, traces obtained by AVO processing. 
The rock properties model includes properties of a seal rock and a reservoir rock. The 
model may be based on two half spaces in contact or may consist of a reservoir rock 
15 sandwiched between two half spaces of seal rock. In the latter case, a tuning curve is 
determined based on an overburden model. 

[0018] Determination of the PDF for the fluid may be based on a convolutional model. 
There is wide latitude in defining the wavelets for the convolutional model. These 
20 include a wavelet derived from a bandpass filter, a Berlage wavelet, (iii) a wavelet 

derived from a Butterworth filter, (iv) a Gabor wavelet, (v) a Gaussian wavelet, (vi) an 
Ormsby wavelet, (vii) a Rayleigh wavelet, and, (viii) a Ricker wavelet. The rock 
properties for the model may be based on trend curves. Compressional wave velocities 
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for the models are directly derived from trend curves, whereas shear velocities and 
densities are based on deviations from expected values based on correlations with 
compressional wave velocities. 

5 [0020] The fluid property may be one or more of a fluid modulus, a density, and, a fluid 
saturation. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[0020] The file of this patent contains at least one drawing executed in color. Copies of 
10 this patent with color drawings(s) will be provided by the Patent and Trademark Office 
upon request and payment of the necessary fee. The present invention will be better 
understood by referring to the following detailed description and the attached drawings in 
which: 

FIG. 1 is a flow chart showing an overview of the preferred embodiment of the invention; 
15 FIG. 2 is a schematic illustration showing the marginal probability density functions for 
an attribute for the data and three different models. 

FIG. 3 is an example of seismic data from the Gulf of Mexico including a horizon of 
interest.. 

FIG. 4 shows a seismic structure map for the horizon of interest. 
20 FIG. 5 is shows the seismic amplitude for the horizon of interest for the area shown in 
Fig. 4. 

FIG. 6 is an example of a PDF for the fluid bulk modulus for a gas filled interval and a 
brine filled interval 
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FIG. 7 is a plot of the most likely fluid modulus over the area shown in Fig. 4. 

FIG. 8 is an example of a PDF for the density for a gas filled interval and a brine filled 

interval. 

5 DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

[0020] Seismic reflection data provides information about reflection coefficients. 
Depending on source and receiver types, one can obtain compressional, shear or 
converted wave reflection coefficients. Reflection coefficients also vary with the 
incidence angle of the reflections. These various reflection coefficients and their angle 
10 dependence provide a variety of information about the two rocks at whose contact the 
reflection occurs. 

[0021] But reflection coefficients have one big limitation: they contain information about 
the pore fluids only in an indirect fashion. It is the pore fluid properties (modulus, 
15 density and hydrocarbon saturation) that we would like to have, because they are directly 
related to the economics of the prospect.. 

[0022] Unfortunately, we know that deriving pore fluid properties from seismic data is 
usually an ill-conditioned process. A variety of rock physics parameters must be known 
20 to go from reflection coefficients to pore fluid properties, and these are not known with 
great precision. In fact, variations in these parameters can mimic the effects of pore fluid 
variations on the reflection coefficients. 
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[0023] The present invention is a method that estimates the pore fluid parameters, and, at 
the same time, also quantifies the uncertainty in the estimates, given our uncertainties in 
the rock properties and the seismic measurements. This is done by producing a 
probability density function (PDF) which gives the probability that the observed seismic 
5 is consistent with a given pore fluid property. 

[0024] This PDF contains all the information about the pore fluid property that can be 
obtained from the seismic measurement, given the uncertainties in our knowledge. The 
PDF contains the most likely value of the pore fluid property. The PDF also quantifies 
10 the uncertainty in the best estimate. The PDF can also be integrated to give the 
probability for a range of values (e.g., those that are likely to reflect pay). 

[0025] The PDF also quantifies how much information is gained from the quantitative 
seismic information. This process can incorporate all the rock properties, log data, etc. 
15 that we have about the problem. If new information is added, we can incorporate it to 
obtain a new (sharper) PDF. Thus we can quantify the importance of new information, 
be it rock properties data or additional seismic measurements. 

[0026] Turning to Fig. 1, there are two broad streams of processing. The left hand side 
20 depicts processing that is used for generation of synthetic data. This may be referred to 
as the simulation. The right hand side depicts the processing done to real data. The 
simulation process starts with a definition of rock properties 101. These are 
characterized statistically and are discussed in detail later in the application. Each 
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parameter is described by a statistical distribution. A plurality of realizations (typically, 
tens of thousands) are drawn from these distributions 103. The variables are chosen to be 
(nearly) uncorrelated. Associated with each of these realizations is a value of a pore fluid 
property chosen to be at the center of a bin. The pore fluid property may be one or more 
5 of the fluid modulus, the fluid density, the fluid saturation (for multi-component fluids), 
etc. For each realization a value of compressional and shear velocities and density (Vp, 
Vs and p) are determined for the caprock (seal) and reservoir rock. This may be done 
for both a test region (where fluid properties are to be determined), and a reference 
region where fluid properties are known (in a statistical sense). 

10 

[0027] For the realizations, equidistributed sequences are used over specified ranges. 
This makes it possible to use fewer realizations than would be necessary if random 
number sequences were used for defining the rock properties. Using the values of V p , 
V s , and p, reflection coefficients are determined. Mathematically, they are obtained by a 

15 solution of Zoeppritz's equations, (given in a classic paper in 1919). A convenient form 
of the equations is given in The Rock Physics Handbook and is reproduced in the 
Appendix. In a preferred embodiment of the invention, the isotropic form of Zoeppritz's 
equations. This is not intended to be a limitation, and the method of the present 
invention could also be used with an anisotropic form of Zoeppritz's equations given, for 

20 example, in The Rock Physics Handbook. 



[0028] In one embodiment of the invention (called the impedance model), reservoir 
impedances are determined. In this embodiment, the seismic data must be inverted to 
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determine reservoir impedance. In a second embodiment of the invention, an assumption 
is made that there are two half spaces in contact, one a seal, the other a reservoir. The 
selected horizon may be either the top or the bottom of the reservoir. This may be 
referred to as the two half-spaces model. The corresponding seismic data must also 
5 satisfy this assumption. In a third embodiment of the invention, a reservoir layer is 
interposed between two half spaces of caprock, whose properties are assumed to be the 
same above and below the reservoir. This may be referred to as a sandwich model. The 
two half-spaces model and the sandwich model both assumes a convolutional model. 
For the two half-spaces model, a convolutional model may be used. However, in a 

10 preferred embodiment of the invention As would be known to those versed in the art, in 
a convolutional model, the seismic response is modeled by convolving a wavelet with a 
reflectivity sequence. For the sandwich model, a tuning curve (derived from the source 
wavelet) is used to replace the two reflections (one from the top and one from the bottom 
of the reservoir) with an equivalent single reflection, and the rest of the processing 

15 proceeds as for the two half-spaces model. 



[0029] In one embodiment of the invention (called the impedance model), reservoir 
impedances are determined. In this embodiment, the seismic data is inverted to 
determine reservoir impedance. In a second embodiment of the invention, an assumption 
20 is made that there are two half spaces in contact, one a seal, the other a reservoir. The 
selected horizon may be either the top or the bottom of the reservoir. This may be 
referred to as the two half-spaces model. The corresponding seismic data must also 
satisfy this assumption. In a third embodiment of the invention, a reservoir layer is 
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interposed between two half spaces of caprock, whose properties are assumed to be the 
same above and below the reservoir. This may be referred to as a sandwich model. The 
two half-spaces model and the sandwich model both assume a convolutional model. As 
would be known to those versed in the art, in a convolutional model, the seismic response 
5 is modeled by convolving a wavelet with a reflectivity sequence. For the two half-space 
model, the reflection amplitude is assumed to be a reflection coefficient multiplied by the 
amplitude of the wavelet at the time picked (e.g. at the peak of the wavelet). It is 
assumed on input that either the wavelet amplitude scale factor is removed or that the 
input will be normalized internally. For the sandwich model, a tuning curve (derived 
10 from the source wavelet) is used to replace the two reflections (one from the top and one 
from the bottom of the reservoir) with an equivalent single reflection, and the rest of the 
processing proceeds as for the two half-spaces model. 



[0029] Returning now to Fig. 1, using the various realizations of the formation plus fluid 
15 models, seismic attributes are determined 105. The seismic attributes are based on the 
impedance (for the impedance model) and on reflection coefficients for the two half- 
spaces model and the sandwich model. The synthetic data are then normalized 107. 



[0030] The processing of the real seismic data (on the right side of Fig. 1) starts with the 
20 seismic data 109. It is input as the values on a mapped horizon. Next, the errors in the 
data are estimated 111. It is assumed that the errors are distributed normally, with the 
mean value(s) being the input value(s). The standard deviation(s) may be input by the 
user and are a specified parameter. In one embodiment of the invention, the PDFs are 
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assumed to be laterally invariant. In an alternate embodiment of the invention, the PDFs 
could vary laterally. Alternatively, the data errors may be estimated data using pre-stack 
data. The data are then normalized 113 using the reference region, if the absolute scale 
of the seismic data is not available. Use of the reference region is as in prior art, i.e., the 
5 values of the attributes on real seismic data in the reference region (where the rock and 
fluid properties are known) are used to normalize the seismic data elsewhere. 

[0031] The (normalized) synthetic and (normalized) real data are then compared 115. 
Based on this, a probability that the real data results from a given value of the fluid 
10 properties is determined. The determination of this probability is depicted schematically 
in Fig. 2. 

[0032] Turning now to Fig. 2, the curve 205 depicts the PDF of a selected attribute for a 
single point on a seismic horizon. The center of the distribution is the measured value of 

15 attribute on the real data at the selected point. Also shown in Fig. 3 are different PDFs 
(201, 202, 203) on the synthetic data for a different selected value of the pore fluid 
property that is to be determined. Three curves are depicted purely for illustrative 
purposes and is not intended to be a limitation. The hatched area 215 is a measure of the 
probability that the observed value corresponds to the fluid property that was used to 

20 generate the curve 202. Similar hatched areas can be determined for the other 

hypothetical values of the fluid properties that gave the curves 201 and 203. For the 
example shown, the fluid property values that gave rise to the curve 202 is the most 
likely value for the real data. The area under both curves (hatched region) is a measure 
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of the probability that the observed value is the result of the chosen pore fluid property. 
Using this, an a posteriori probability for the fluid property (e.g., fluid modulus, the fluid 
density, the fluid saturation) in question can be calculated 

5 [0033] The example given above in Fig, 2 dealt with a PDF for a single attribute. This 
methodology can be extended to the measurement of multiple attributes, in which case, a 
joint PDF would be used. 

[0044] One of the important steps in the invention is the realization of the rock properties 
10 model and the effects of fluid modulii (103 and 105) in Fig. 1. Some important 
considerations in these two subprocesses are discussed next. 

[0045] One of the principles used in the modeling is that the P and S velocities of rocks 
should trend between the velocities of the mineral grains in the limit of low porosity and 

15 the values for a mineral -pore fluid suspension in the limit of high porosity. See Nur et al 
(1995). This idea is based on the observation that for most porous materials, there is a 
critical porosity, <p c , that separates their mechanical and acoustic behavior into two 
distinct domains. For porosities lower than <p c the mineral grains are load bearing, 
whereas for porosities greater than <p c the rock falls apart and becomes a suspension in 

20 which the fluid phase is load bearing. Based on this, for (f> > <p c , the effective bulk 

modulii and the shear modulii can be estimated quite accurately using the Reuss average: 

K- R l =(l-(/>)K; l + <t>K- fl l ; ju R =0 (1) 
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where K R . K Qy and K fl are the bulk modulii of the rock, matrix and fluid respectively, 
and fi R is the shear modulus of the rock. 



[0046] In load bearing rocks, the modulii decrease rapidly from the mineral values at 
zero porosity to the suspension values at the critical porosity. Nur found that this 
dependence can often be approximated with a straight line when expressed p V 2 versus 
porosity. For dry rocks, the bulk and shear modulii can then be expressed as: 



K dry ~ K gr 



t^dry ftgr 



Yc J 



(2) 



In a preferred embodiment of the invention, a modified form of the equations given by 
10 Nur are used: 



K, , — K n 

dry gr 



M dry =M, 



\ TC J 



(2a) 



where r is an exponent. It has been found that this model gives a better result than the 



GTI-1539 



17 



one given by Nur. 

[0047] In a preferred embodiment of the invention, the rock properties are specified as 
follows. 

SEAL ROCK PROPERTIES are determined from user specified trend curves from the 
5 relations: 

V s = Trend function (V p ) 

p = Trend function (V p ) (3). 

Such trend curves are well known in the art and are usually compiled from measurements 
made in wells on a local scale, a prospect scale, or on a basin-wide scale. Usually, the 

10 sampling of the subsurface by wells is rather sparse on a local scale; hence trend curves 
are preferably obtained on a prospect scale or a basin-wide scale. In a preferred 
embodiment of the invention, a trend function of V is specified. The plurality of 
realizations are based on perturbations of V p about this trend curve. As is well known, 
there is a strong correlation between V s and V p . Hence in a preferred embodiment of the 

15 invention, the plurality of realizations are specified in terms of perturbations of V s about 
the trend. Thus, if S V p is a perturbation of V p and if A V 5 is the expected deviation in V s 
for a perturbation 3 V p in V p , then a perturbation 6 V s is specified as additional to AV S . 
This ensures that the realizations comprise uncorrected perturbations. A similar step is 
taken with respect to density, which is also known to be highly correlated with V p . 

20 

RESERVOIR ROCK PROPERTIES are determined using the well known Gassman 
equations and trend curves as: 

K dry = K grain ( 1 " <P / <Pc) r ( 4 ) (critical porosity model) 
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p 

M 

\IK, 



■fluid 



r brine 

Pgas 

Pfluid 



= PgrainW " <P) + Pfluid <P 

= trend function (K dry ) 
= (1 - S g )IK brine + S g /K gas 
= (l-S 0 )/K brine + S 0 /K oil 
= Trend, function (K brine ) 
= Trend.function(K^) 
= (1 - S g )p brim + S g p g 

= (1 -S 0 )Pbrine + SoPg 



(5) 

(6) (for gas saturated fluid 
(for oil saturated fluid) 

(7) 
(8) 

(9) (gas saturated fluid) 
(oil saturated fluid) 



K = K^. + 



(! _ K dry I K grain ) 



dry 



{{^-(l>-K dry IK grain )lK grain +<fi/K Jluid ) 



(10) (from Gassman) 



10 V p =A(* + 4u73)/p) (11) 

V s = v/(uVp) (12) 



The single phase fluid density p from the K trend function is: 
p = .375 K 1/2 for K <. 1 (Gas region) and 

15 p =J5K m forK>,15 (13) (Liquid region). 

These equations are based on the density having units of g/cm 3 and the modulii being in 
GPa. Expressions for other units are easy to derive. For the reservoir the stochastic 
variables are K^, p^,,, porosity, critical porosity, exponent r (see above), change of 
shear modulus from trend, brine modulus, brine density deviation from trend. 

20 
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[0048] Another important consideration in Fig. 1 is the normalization of the data 107, 
113. In a preferred embodiment of the invention, a normalization of the form: 

/y(D test 2 +D ref 2 ). (14) 
is used. This choice of form prevents the large values when either data value is near 
zero. The reference region will normally contain many points. A small but 
representative subset of these is chosen and processed for each case and the resulting 
PDFs. are averaged. This is the equation for generating the normalized input data, D n01 
from the input data value in the test region, D test , and the input value in the reference 
region, D refi 



[0049] The invention can use data (and attributes) from a variety of formats. Some of the 
more common options for the data are: 

1. Both compressional and shear wave data may be used. 

2. The data may be stacked or unstacked. 

15 3. Mode converted data (PS or SP) are acceptable. 

4. Data may have been processed using commonly used AVO analysis. 
Specifically, A, B, C, D and E traces may be used as commonly defined. 

5. Data may have been processed to give reflectivities R p> R s and impedances Z p and 

20 6. Broad flexibility exists in the wavelet definition. Among the commonly used 
wavelets are (i) a wavelet derived from a bandpass filter, (ii) Berlage, (iii) a 
wavelet derived from a Butterworth filter, (iv) a Gabor wavelet, (v) a Gaussian 
wavelet, (vi) an Ormsby wavelet, (vii) a Rayleigh wavelet, and, (viii) a Ricker 
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wavelet. These wavelets would be known to those versed in the art and are not 
discussed further. 

7. When prestack data are analyzed, an overburden model may be defined in terms 
of velocity and density to make appropriate calculations of the angle of incidence. 
5 These are only a sampling of some of the more important parameters that may be 
defined. 

[0050] With the above in mind, the attributes used in the analysis may be selected from 
the list below. Reference is made to prior art methods of analysis of offset dependent 
10 seismic traces wherein: 

R(0) = A + B sin 2 (6) + C sin 2 (0) tan 2 (0) (15), and 

R(0) = D sin(0) + E sin 3 (0) (16), 

where R(0) is a seismic trace corresponding to an angle of incidence 0, and A, B, C, D 
and E are fitting parameters. As would be known to those versed in the art, eq. 15 which 
15 is a Taylor series expansion in even powers of Sin0 would be used for either P-P 

(incident compressional wave, reflected compressional wave) or for S-S (incident shear 
wave, reflected shear wave) data, while eq, (16) which is a Taylor series expansion in 
odd powers of Sin0 would be used for (P-S) or (S-P) data. An example of such a method 
is given (for P-P data) in US Patent 4,995,007 to Corcoran et al. 

20 

The measured attributes are then one or more of: 

1. Z p the compressional wave impedance, 

2. Z s the shear wave impedance (shear), 
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3. P-Wave stack amplitude, 

4. S-Wave stack amplitude 

5. PS converted stack amplitude, 

6. SP converted stack amplitude, 
5 7. AVO attribute A (P-wave), 

8. AVO attribute D (PS Converted), 

9. AVO attribute A (S-Wave), 

10. AVO attribute D (SP Converted), 

11. Reflectivity (P-P), 
10 12. Reflectivity (S-S), 

13. Impedance (P- and S-), 

14. Reflectivity (P- andS-), 

15. Stack amplitude (P- and S-), 

16. AVO attributes A and B (P- waves), 

15 17. AVO attributes D and E (PS converted waves), 

18. AVO attributes A and B (S- wave), 

19. AVO attributes D and E (SP converted waves), 

20. AVO attributes A, B and C (P-wave) 

21. AVO attributes A, B and C (S-Waves) 

20 22. Reflectivities (P and S), fractional change in density, and 
average (Vp/Vs) 2 
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[0051] Turning next to Fig. 3, an example of a seismic line is shown. This is a line of P- 
P data. Indicated by the arrow is a horizon of interest for which an estimate of the fluid 
content is to be determined. Fig. 4 shows a seismic structure map for an area of interest. 
Contours of equal reflection times of the horizon of interest are indicated. There are 
5 several faults in the area, indicated by the dark blue lines. Attention is drawn to the fault 
blocks indicated by 301, 302, 303 and 305 which will be discussed further below. In 
addition to the contours, the color is an indication of the depth of the horizon (the scale 
ranges from 2.35 to 2.65 seconds). Several wells have been drilled in the area as shown 
by the red dots in Fig. 4. Data from these wells may be used to define the rock 
10 properties and fluid properties at known calibration points in the area. 

[0052] Fig. 5 is a display of the seismic amplitude over the region of interest. As can be 
seen, the amplitudes are the highest in block 301 as indicated by the red color, while they 
are the lowest in block 305 ad indicated by the blue color. Intermediate values of the 
15 amplitude are noted in blocks 303 and 303. Qualitatively, the high amplitudes (red 

values in Fig. 5 and the strong amplitudes seen in Fig. 3 at the horizon of interest) are 
indicative of gas. The rest of the analysis used only the amplitude of the P-P stack as an 
attribute. 

20 [0053] Turning to Fig. 6, a plot is shown of the PDF that would results from a gas filled 
interval and a brine filled interval. Not surprisingly, a brine filled interval has a high 
PDF at a large value of the fluid bulk modulus while a gas filled zone has a high PDF at a 
small value of modulus. Fig. 7 shows the most likely value of the fluid modulus (in 
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GPa) over the region of interest. As can be seen, the crest of the structure in block 301 
has the lowest likely value fluid modulus (indicative of gas) while the block 305 has the 
highest likely value of fluid modulus (indicative of brine); blocks 302 and 303 have 
intermediate values. One might infer that the structure has gas. However, any 
5 conclusion based on the data that the reservoir has a high gas saturation would be 
unjustified for reasons discussed next. 

[0054] As would be known to those versed in the art, eq. (6) shows that low fluid 
modulus is not particularly diagnostic of high gas saturation. There is a common 

10 occurrence of what is called low saturation gas (LSG) wherein strong seismic reflections 
and low fluid modulus result from low values of S g . It is important to get an estimate of 
the density as well. Eq. (9) shows that the fluid density is sensitive to the gas saturation. 
However, in the present example, just using the seismic amplitude alone provides very 
little information about the fluid density. This can be seen from Fig. 8 which shows a 

15 plot of the PDF for gas and brine filled intervals that could be obtained using the 

amplitude as the only attribute. Fig. 8 shows that no inference can be drawn about the 
fluid density (or gas saturation) from the seismic data when the P-P amplitude is the only 
seismic attributes. Other attributes must be included. 

20 [0055] In a preferred embodiment of the invention, P-P seismic data are used, the 

attributes measured are A and B for the data, and the fluid parameters for which the PDFs 
are determined are the fluid modulus and the fluid density. 
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[0055] The foregoing description has been limited to specific embodiments of this 
invention. It will be apparent, however, that variations and modifications may be made 
to the disclosed embodiments, with the attainment of some or all of the advantages of the 
invention. Therefore, it is the object of the appended claims to cover all such variations 
5 and modifications as come within the true spirit and scope of the invention. 



APPENDIX 

The general formulation for a plane wave incident on an interface between two half 
spaces is given by 



pIpI slpl 
pisl sisi 
pIpI si pI 
pisi sisi 



pi pi si pV 
pisl sisi 

Pi Pi Si Pi 

pisi sisi 



= M~ l N 



(A-l) 



where the first arrow defines the direction of the incident wave and the second arrow 

defines the direction of the reflected wave and where the matrices M and N are given by 

-sintf, -cos$ sin# 2 cos^ 

cos#, -sin$ cos6 2 -sin<p 2 

M ~ 2p x V sX sin <j\ cos 0, p, V st (1 - 2 sin 2 $ ) 2p 2 V s2 sin 0cos 0 p 2 V s2 (1 - 2 sin 2 fa ) 
_-/V pI (l-2sin 2 #) P 1 V J1 sin2^ p 2 V p2 (l-2sin 2 0) - p 2 V l2 sin 2^ 
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and 



N = 



sin#j 
cosfl 



cos$ 
sin$ 



-sin 0 2 
cos ft 



-cos 02 
- sin (f> 2 



2p x V sl sin <j\ cos 0 X p,V J1 (l-2sin 2 </\) 2p 2 V s2 sin0cos0 p 2 V s2 (l-2sin 2 ^) 
flV pI (1-2 sin 2 ft) -fll/ fI sin 2$ (1-2 sin 2 <p) p 2 V s2 sin20 2 



(A-2) 
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In the equations, the angles 0 and <j) refer to the angles of P- and S- waves respectively in 
a first and second medium separated by a plane interface, V p , V s and p refer to P- and S- 
15 velocities and densities and the subscripts refer to the first and second medium. Various 
approximations of solutions of Zoeppritz's equations have been made. 



For normal incidence, the solutions are of the form: 



20 For small angles of incidence of an incident P- wave, the reflection coefficient R pp ( 6), of 
the reflected P- wave i.e, the Pi PI term above, has a value that is approximately given 
by 
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R.J0) *A + Bsin 2 d 



(A- 4) 



where the coefficient B depends upon the difference in shear velocity across the 
reflecting interface. Similarly, the reflection coefficient R ps ( 6) of the reflected shear 
5 wave for an incident P- wave is approximately given by a relation of the type 

R ps (6) = Csin(6) (A-5) 

Qualitatively, the amplitude of the reflected P- wave depends upon the angle of incidence 
10 of the P- waves at the interface. 
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